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We investigate an application of twisted boundary conditions for the study of low-energy hadron- 
hadron interactions with Liischer's finite size method. This allows us to calculate the phase shifts 
for elastic scattering of two hadrons at any small value of the scattering momentum even in a finite 
volume. We then can extract model-independent information of low-energy scattering parameters, 
such as the scattering length, the effective range, and the effective volume from, the 5-wave and 
P-wave scattering phase shifts through the effective range expansion. This approach also enables 
us to examine the existence of near-threshold and narrow resonance states, whose characteristics 
are observed in many of newly discovered charmonium-like XYZ mesons. As a simple example, we 
demonstrate our method for low-energy J/ip-(j> scatterings to search for Y(4140) resonance using 
2+1 flavor PACS-CS gauge configurations at the lightest pion mass, m w — 156 MeV. 



I. INTRODUCTION 

In past several years, properties of hadronic interac- 
tions have been extensively studied in lattice QCD sim- 
ulations based on Liischer's finite size method, which 
is proposed as a general method for computing low- 
energy scattering phase shifts of two particles in finite 
volume [TJ [2] . Various meson-meson, meson-baryon and 
baryon-baryon scattering lengths, which are relevant for 
describing low-energy scattering processes, have been 
successfully calculated in lattice QCD within this ap- 
proach [3]. 

The original Liischer method, which is considered for 
zero total momentum P = of the two-particle system 
in a symmetric box of size I? with periodic boundary 
conditions, has been developed in various ways [4] . In 
particular, there are several extensions of its formula for 
computing the scattering phase shifts at more values of 
lower scattering momenta in a given volume. From a 
practical viewpoint, accessible values of the phase shift 
on the lattice are limited due to discrete momenta, ap- 
proximately, in units of 2ir/L. To increase accessible mo- 
menta in a single volume, the formalism was extended 
to moving frames, where the scattering particles have 
nonzero total momentum [7HH], and also general- 

ized in an asymmetric box of size {r)iL) x (772-EO x L, where 
the degeneracy of low-lying modes in three-dimensional 
momentum space can be resolved for 771 , 772 ^ 1 [101 lllj . 
Alternatively, we notice that an idea of twisted boundary 
conditions is quite useful for studying the hadron-hadron 
interaction at low energies through Liischer's finite size 
method as originally proposed by Bedaque |12) . 

Since the twisted boundary conditions allow us to eval- 
uate scattering phase shifts of the two-particle system 



* sho@rcnp.osaka-u.ac.jp 
tssasaki@nucl.phys.tohoku.ac.jp 



at any small value of the scattering momentum even 
in a finite box, detailed information of the low-energy 
interaction, which is represented by the lower partial- 
wave phase shifts at low energies, can be easily obtained 
through an extended Liischer formula for twisted bound- 
ary conditions. In general, the quantity of k 2l+1 cot 6i(k), 
where Si(k) and k denote the phase shift of the Ith partial 
wave and the absolute value of the scattering momentum 
(k = \k\), can be expanded in a power series of the scat- 
tering momentum squared in the vicinity of the threshold 
as 

/fc 2/+1 cot(5 ; (fc) = - + 1 r l k 2 + C(fc 4 ), (1) 
ai 2 

which is called the effective-range expansion [13] . Model- 
independent information of the low-energy interaction is 
encoded in a small set of parameters, e.g., scattering 
length ao and effective range tq for an S wave (I = 0). 
Therefore, we can determine such scattering parameters 
through the k 2 dependence of the scattering phase shifts 
near the threshold using the novel trick of twisted bound- 
ary conditions |14| . 

Furthermore, we are able to apply this method to in- 
vestigate newly observed narrow resonances in the heavy 
sector. Recently, many charmonium- and bottomonium- 
like resonances, the so-called XYZ [HI [16] and Y^Zi, 
[TTH20] resonances, are reported in several experiments. 
Among them, some resonances are observed near two- 
particle thresholds, and their widths are quite narrow 
as compared to typical hadron resonances. As a sim- 
ple application for the new approach, we study a low- 
energy scattering of two mesons and search for a nar- 
row resonance near the threshold. The J/tp-cj) channel 
is considered to be an appropriate research target, since 
three narrow resonances have been reported in recent ex- 
periments, namely, F(4140) and F(4274) by the CDF 
Collaboration [2]] [22] and AT (4350) by the Belle Collab- 
oration [23]. Interestingly, these resonances seem to be 
relatively stable despite being above open charm thresh- 
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olds, since the upper bound of their widths is less than 
10-30 MeV. In particular, F(4140) is located close to 
the J/4>-4> threshold. The F(4140) resonance was first 
reported by the CDF Collaboration with 3.8a statistics 
in 2009 [H] . The signal is observed in the invariant mass 
of the J I pairs of the decay B + — > J / ip<fiK + from pp 
collisions at y/s = 1.96 TeV. A preliminary update of the 
CDF Collaboration analysis leads to the observation of 
the Y"(4140) with a statistical significance of more than 
5a [22]. The mass and width are 4143.4^|q ± 0.6 MeV 
and 15.3l6° 1 4 ±2.5 MeV, respectively [22]. Although the 
observed mass is much higher than the DD threshold, 
the V(4140) resonance has a very narrow width. This 
observation suggests that the V(4140) scarcely couples to 
open charm channels such as DD and also D+D~. On 
the other hand, the Belle and LHCb collaborations have 
not yet found the Y(4140) in their experiments [2"31 124) . 
Although the Y"(4140) might have very interesting prop- 
erties, its existence is still controversial experimentally. 
Our analysis of low-energy J/ip-tfr scatterings under the 
twisted boundary conditions could give a new insight into 
the Y"(4140) resonance from first principles QCD. 

The paper is organized as follows. In Sec. II, after a 
brief introduction of the twisted boundary condition, we 
show the finite size formula with such particular bound- 
ary conditions. Next, we apply our formula to the low- 
energy J/iJj-(j) scattering in Sec. Ill, and then show our 
results in Sec. IV. Finally, we summarize our study. 

II. THEORETICAL FRAMEWORK 

A. Boundary conditions 

Let us recall the ordinary periodic boundary condition 
in the spatial directions: 

+Le u t) = ^(x,t) (2) 

with the Cartesian unit vector ei along the i axis (i = 
x, y, z), which provides the well-known quantization con- 
dition on three-momentum for the noninteracting case: 
p = {2ir/L)n with a vector of integers n G Z 3 . A typical 
size of the smallest nonzero momentum under periodic 
boundary conditions, e.g., \p m i n \ ~ 2-7r/L ~ 420 MeV for 
L ~ 3 fm and 250 MeV for L ~ 5 fm, is still too large to 
investigate the hadron-hadron scattering at low energies. 

A novel idea, the twisted boundary condition, was pro- 
posed by Bedaque to circumvent this issue [12]. The 



twisted boundary condition is a sort of generalization of 
the periodic boundary condition in the following way: 

^e{x + Lei,t) = e i6i ^ e {x,t), (3) 

where the angle variable 9 is called the twist angle and 
tyg stands for either the elementary or composite field 
operator on which twisted boundary conditions are im- 
posed. Here 9i = corresponds to the ordinary periodic 
boundary condition as described above, while 9i — n cor- 
responds to the antiperiodic boundary condition. Under 
twisted boundary conditions, the discretized momentum 
on the lattice should be modified as 

p= T\ n+ Y,) (4) 

with a twist angle vector 9 = (9 x ,9 y ,9 z ) for the free 
case. In principle, we can have any value of momentum 
on the lattice through the variation of the twist angle, 
continuously. In particular, in this paper, we want to 
emphasize that the lowest Fourier mode, n — (0,0,0), 
can still receive nonzero momentum, which can be set to 
an arbitrary small value even in a fixed spatial extent L, 
using this novel trick. 

Of particular interest is the case of partially twisted 
boundary conditions, where we impose twisted bound- 
ary conditions on the valence quark fields of specific fla- 
vor. Needless to say, if the twisted boundary conditions 
are imposed on the sea quark fields, huge computational 
cost is required to generate a new gauge ensemble for 
each twisted angle. In addition, the authors of Ref. [2"5] 
show that the finite volume correction due to the partially 
twisted conditions is exponentially suppressed as the spa- 
tial extent L increases. Therefore, there is the practical 
advantage of the partially twisted boundary conditions. 

For convenience of explanation, we introduce new fields 
q' , which are transformed from the original fields qg where 
the twisted boundary conditions are imposed as 

q'(x,t) = e- i ^/ L q e (x,t). (5) 

These fields now turn out to satisfy the usual boundary 
conditions as q' (x + Lei) — q'(%)- This redefinition of the 
quark fields affects only the hopping terms that appear 
in lattice fermion actions. For Wilson-type fermions, the 
hopping terms are transformed as 



1e(x) (1 - iu)U x ^8 x+ll<v + (1 + IvdUt-^Sx-w qe(y) 



= ]r q'(x) Ui ~ jMje)s x+M + (i + ^w'l^m 



x-fj,,y 



q'(y) 



(6) 
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FIG. 1: The grid points in momentum space specified by a 
vector r are plotted in the xz plane. For the case of a twisted 
angle vector 6 — (0, 0, 9), the grid points are shifted from the 
origin in the z direction by d z — 9/2-n. The panels show the 
cases of (a) 8 — 0, (b) < 6 < -n and (c) 9 = ir. A center of 
symmetry is clearly lost in the middle panel. However, when 
9 — tt (right panel), the origin is shifted at the midpoint of 
the grid interval, where the inversion center gets positioned. 



with U' x Jfi) 



e ,^a/L Ux ^ and ^ = (q ^ Therefore, 
we can easily calculate the quark propagator subject to 
twisted boundary conditions through the simple replace- 
ment of the link variables {U x ^} by {U'(a}} in the 
hopping terms. 

Using this technique, we can easily construct a 
hadronic interpolating operator for states with finite mo- 
menta. We now consider the traditional meson interpo- 
lating operator as a local bilinear operator, Or(x,t) = 
q'f(x,t)Tq'f,(x,t), where /, /' denote flavor indices with 
Dirac's gamma matrices T, as a simple example. A sim- 
ple summation over the spatial sites on this operator can 
be interpreted as the Fourier transformation to a momen- 
tum representation: 

X> r (f,t) = ^^/(^^rge,/'^*)^'" 5 *' '^ 

= O r {p,t), (7) 

wherep= (9 qf —0q,)/L. Unless the same twisted bound- 
ary conditions are imposed on both quark and antiquark 
operators, the resulting meson operator does receive fi- 
nite momentum because 9 qf ^ 6q , . This technique is 
widely used for flavorful mesons (/ ^ f) in several lat- 
tice analyses 



B. Finite size formula 

Let us consider the finite size formula under partially 
twisted boundary conditions. In this study, we only con- 
sider the center-of-mass (CM) system with zero total mo- 
mentum P = 0, so that the CM system and the labo- 
ratory system coincide. Therefore, we do not have to 
consider the Lorentz boost from the laboratory frame to 
the CM frame in our approach, unlike recent works that 
considered moving frames [25H3"3"] . 

First of all, the Liischer finite size formula provides 
a relation between finite volume corrections to the two- 
particle spectrum and the physical scattering phase shift 
through a consideration of the relative two-particle wave 



function ^(x,t) in the CM frame. Since ^(x, t) defined 
in the CM frame should be independent of the relative 
time t, we simply omit the argument for the relative time 
hereafter. Here the wave function ^(x) is supposed to 
satisfy the Hclmholtz equation in an outer range, where 
an interaction potential of finite range vanishes [2]. The 
Helmholtz equation is obtained with a set of wave num- 
bers associated with discrete relative momenta fc 2 in a 
finite box L 3 : 



(A + k 2 )^(x) = 0, 



(8) 



where fc 2 corresponds to the scattering momentum. So- 
lutions of the Helmholtz equation are supposed to satisfy 
the twisted boundary conditions ([3]). Here we introduce 
the Green function, which obeys the twisted boundary 
conditions: 



G e {x,k 2 ) 



" 3 E 

per,- ' 



p 2 



fc 2 ' 



(9) 



where the momenta p are the elements of = {p\p — 

~H + ~,n e Z 3 }. This function is of course a solu- 
tion of the Helmholtz equation for x ^ (mod L). Fur- 
ther solutions may be obtained as its derivatives by using 
the harmonic polynomials, yi m {x) = x l Yi m (9,(f>) defined 
with spherical coordinates, x — (x, 9, <f)): 



Gf m (x,k 2 ) = y lm (V)G e (x, fc 2 ). 



(10) 



These functions form a complete basis of the singular pe- 
riodic solutions of the Helmholtz equation, which are sup- 
posed to have degree A. The solutions can be represented 
by a linear combination of the functions Gf m (x, fc 2 ) with 
arbitrary coefficients vi m : 



A 

E 

1=0 



E v imG e lm {x,k 2 ), 



(11) 



which satisfies lim^^o |a; A+1 ^E'(x) | < oo. 

Here it is worth mentioning that symmetry proper- 
ties of G (x, fc 2 ) are not the same as the space group of 
the cubic lattice, namely the octahedral group Oh, when 
9 =/= 0. This is simply because the finite momentum is in- 
duced by the twisted boundary conditions, and therefore 
the symmetry of the reciprocal lattice is reduced to the 
subgroup of the full space group. Different partial-wave 
contributions are highly mixed in the expansion expres- 
sion of G lm (x, fc 2 ) in terms of spherical harmonics due to 
further reduction of the symmetry. 

When we take 9 = 0, the cubic symmetry is surely sat- 
isfied in both the space lattice and the reciprocal lattice. 
The partial-wave mixing is maximally suppressed. In 
particular, even-l and odd-Z waves do not mix with each 
other because of a center of inversion symmetry in the cu- 
bic group. The irreducible representation A\ g of the cu- 
bic symmetry contains partial waves of / = 0, 4, 6, 8, • • • < 
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A so that the wave function projected in the irrep Ai g 
may still mix I = and I = 4 waves at low energies [TJ [2] ■ 
As far as we consider the S'-wave phase shift S near the 
threshold, the higher partial-wave (I > 4) contributions 
are safely ignored. One may then choose an angular mo- 
mentum cutoff as A = for the irrep A\. The original 
Liischer finite size formula for the S'-wave phase shift S 
is then given as the determinant of a 1 x 1 matrix [J [3] : 

cotdo(fc) = -^Z 00 (l;q 2 ), (12) 

TT^I z q 

where k and q are the relative momentum of the two- 
particle system in the CM frame and its scaled momen- 
tum defined by q 2 = (^) with the spatial extent L. 
Here the function Z 00 (s; q 2 ) is the generalized zeta func- 
tion, which is formally defined as 

2 -M 2 )=EM (13) 

for s > 3/2, and then an analytic continuation of the 
function is needed from the region s>3/2tos = l[Tjr2]. 

In the case of the twisted boundary condition with 
twist angles < \9\ < tt, however, the cubic symmetry is 
broken into a subgroup symmetry. For instance, when we 
take a twist angle vector 9 = (0,0,9), which is denoted 
symbolically by [001] hereafter, the grid points in mo- 
mentum space are specified by a vector, r = n + d where 
fi e Z 3 and d = (0, 0, 0/2tt). In the case of \9\ ji 0, the 
mesh points are shifted in the z direction by d z — 9/2tt, 
as schematically depicted in Fig[T] Then, the symmetry 
of the reciprocal lattice, whose vector can be given by 
multiplying the vector r by a factor of 2tt/L, is described 
by the little group C^ v . Although the point group C^ v 
does not have a center of inversion symmetry, it is re- 
covered in the case of a special angle, \9\ = tt, where 
the symmetry should be described by the D^h group as 
shown in Fig. [l](c). Similarly, when 9 = (9,9,0) (labeled 
by [110]), the symmetry is further reduced to point group 
C2v (D2h) for < \9\ < tt (\9\ = tt), while the symmetry 
is described by C 3v (D 3d ) for < |0| < tt (\9\ = tt) in the 
case of another twist angle vector 9 = (9,9,9), which is 
denoted symbolically by [111]. 

As mentioned previously, the point groups C nv for 
any finite n do not have a center of inversion symme- 
try. As a consequence, even-Z and odd-Z partial waves 
would be mixed together under twisted boundary condi- 
tions in the range of < \9\ < tt. The possibility of such 
unwanted mixing is simply ignored in early exploratory 
works with twisted boundary conditions (28] [34]. Re- 
cently, the same pathological issue is properly addressed 
in the derivation of Rummukainen-Gottlieb's formula 
(or, equivalently, Luscher's finite size formula in a mov- 
ing frame) for the case of two particles with different 
masses (30H33j , where the symmetry of the system is re- 
duced to the same little groups C nv . 



Twist angle (0,0, 9) (6,0,0) (0,9,9) 

Symmetry dv Ci v C-i V 

Label [001] [110] [111] 

M e S s w oo w 00 woo 

M% P i\^3wio is/6wn i3wio 



M° PP woo + 2w2o woo — W20 — iy/6w22 woo — i2\/6^22 



TABLE I: Matrix elements A4 ab for three different types of 
the twisted angle vector 9. For simplicity, an irrelevant phase 
factor is omitted in the definition of M%p. The explicit ex- 
pression of wi m is given in the text, wn and W22 are complex 
functions, while woo, wio and W20 are all real functions. If 
special properties of Re(mn) = Im(wn) and Re(w22) = are 
taken into account, ton and W22 can be rewritten by a single 
real function as wu = (1 + i)Re(wn) and W22 = ilm(ui22). 

The Rummukainen-Gottlicb finite size formula is 
known as a generalization of the Liischer formula for two- 
particle states calculated in the laboratory frame, where 
the total momentum of two particles is nonzero. The co- 
ordinate in the CM frame can be related to the original 
coordinate in the laboratory frame through the Lorentz 
transformation. In a moving frame approach, the original 
wave function, which is defined in the laboratory frame, 
can be Lorentz boosted into the appropriate CM frame 
with the Lorentz factor 7 [35] to derive the finite size 
formula, which provides a relation between finite volume 
corrections to the two-particle spectrum and the physical 
scattering phase shift. 

In the original work [7], where two degenerate states 
were considered, the factor simply introduces a negative 
sign in the boundary condition for the relative wave func- 
tion. This indicates that the antiperiodic boundary con- 
dition, instead of the periodic boundary condition, is im- 
posed in the direction parallel to a nonvanishing compo- 
nent of the total momentum. However, for the case of two 
particles with different masses, the factor becomes a unit 
complex number, whose argument depends on the mass 
difference [2H1 ED]- Consequently, the boundary condi- 
tion for the relative wave function becomes identical to 
the twisted boundary condition, where the twist angle 
vector is parallel to the total momentum and the size 
of twist angles is associated with the mass difference of 
two-particles. 

For the case of P = 2tt/L(0, 0, 1), Fu has derived the 
generalized Rummukainen- Gottlieb formula for the triv- 
ial Ai sector of C± v as the determinant of a 2 x 2 matrix, 
which is composed of both I = and I = 1 channels [30] . 
Subsequently, the finite size formula for a moving frame 
is further generalized for different irreducible represen- 
tations of the reduced symmetry and also other types 
of nonvanishing total momenta: P = 27r/L(l,l,0) and 
27r/L(l, 1, 1) [3"TH3"3"] . We easily apply this knowledge to 
our twisted boundary cases in the CM frame, where the 
Lorentz boost factor is unity, 7 = 1. 

The finite size formulas for P 7^ kinematics [3T)l - f3"3"] 
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can be easily translated into the Liischer finite size for- 
mula given under the partial twisted boundary condi- 
tions. Following the original expression for a moving 
frame [30H33j , when we take an assumption that <5;>i = 0, 
the finite size formulas for the A\ irrep are given as the 
following determinant 



cot So(k) 

M% P {q) 



M% s {q) M S SP {q)^ 

cot 5 1 (k) - M PP (q) 



0, (14) 



which properly takes into account the mixing of S- and 
P-wave phase shifts in one relation. The ma trix elements 
M s ss , M e SP and M e PP that appeared in Eq.( 14) for three 
types of the twist angle vector are given in Table U where 
we use the shorthand notation 



wim(q) 



7r 3 / 2 y /2l + lq l + 



(15) 



The generalized zeta function Zf m (l;q 2 ) with the twist 
angle vector is defined as 



zL(s;q 2 



rer ff 



y,, 



(r 2 — q 2 ) s 



(16) 



where the lattice grids f in the momentum space are the 
elements of Tg — {r\r = ft + ft <G Z 3 }. For numer- 
ical evaluation of Zf m (s;q 2 ), we use a rapid convergent 
integral expression found in Appendix A of Ref . [36 . 

The finite size formulas for three types of the twist 
angle vector are superficially quite different from each 
other, and this gives rise to some anxiety about the pos- 
sibility that there is no unique solution of the phase shift 
at given k 2 . However, we have numerically verified that 

110]/ A _, , ,,[1111 



when \q\ « 1/2, Ml b 



whereas when 1/2 < \q\ < \/2/2, M l ab 



l (q)^M^(q). 



Here we note that when we restrict ourselves to the region 
\9\ < 7T, such a region of the twist angle gives the upper 
bound of the allowed kinematical region as \q\ < \ for the 



case of [001], \q\ 



< 



v_2 



for the case of [110], and \q\ 



< 



for the case of [111], respectively. Therefore, three finite 
size formulas could approximately provide the identical 
values of the S-wave and P-wave phase shifts obtained 
at the same scattering momenta within systematic un- 
certainties stemming from the reduction of the lattice 
rotational symmetry from the cubic symmetry Oh to dif- 
ferent little groups C nv . 

As mentioned previously, these modified finite size for- 
mulas show the mixing of S and P waves. However, at 
\6\ = 7T, the parity symmetry is restored, and the 5-wave 
and P-wave phase shifts are disentangled owing to the 
group theoretical constraint of A4g P = 0. In this case, 
Eq. ( 14 ) gives a similar finite size formula to the original 



Rummukainen-Gottlieb one [7] without the Lorentz fac- 
tor (since 7=1). On the other hand, when 9 = (0, 0, 0), 
Eq. ( 14 ) properly reduces the original Liischer finite size 



formula ( 12 ) due to the fact that -M e SP — 0. Here we ne- 
glect higher partial-wave contributions above the D wave 
(I > 2). Such contributions should be kinematically sup- 
pressed as far as we consider low-energy hadron-hadron 
scattering near the threshold. 



III. SIMULATION DETAILS 
A. Lattice setup 

We apply the Liischer finite size method that is gener- 
alized under the partially twisted boundary conditions to 
explore J/ip-ip scattering at low energies. A narrow reso- 
nance Y(4140) that appeared in a J/ip-<fr decay mode has 
been reported by the CDF Collaboration. If the F(4140) 
state really exists, we would directly observe a shape res- 
onance in the J/ip-<fi scattering using lattice QCD at the 
physical point. We thus have performed dynamical lat- 
tice QCD simulations on a lattice L 3 x T = 32 3 x 64 
with 2+1 flavor PACS-CS gauge configurations, where 
the simulated pion mass is closest to the physical point 
as — 156(7) MeV 37 . Simulation parameters of 
PACS-CS gauge configurations are summarized in Ta- 
ble HU Our results are analyzed on all 198 gauge config- 
urations, which are available through the International 
Lattice Data Grid and the Japan Lattice Data Grid |38j . 

We use nonperturbatively improved clover fermions 
for the strange quark (s) and a relativistic heavy quark 
(RHQ) action for the charm quark (c). The RHQ action 
is a variant of the Fermilab approach [55], which can 
remove large discretization errors for heavy quarks. Pa- 
rameters of clover fermions and the RHQ action used in 
this work are listed in Table |Hl| Although the simulated 
strange quarks are slightly off the physical point, the pa- 
rameters are chosen to be equal to those of the strange 
sea quark used in gauge field generation. For the charm 
quark, adequate RHQ parameters defined in a Tsukuba- 
type action [301 EI] were calibrated to reproduce the ex- 
perimental spin-averaged mass of a IS" charmonium state 
in Ref. 142]. 



B. Interpolating operators 

In general, the Fourier transform of hadron interpolat- 
ing operators or the summation over all spatial points 
on the operators like Eq. ^ at the source time slice 
is rather expensive from a computational point of view. 
The stochastic techniques are often used for this purpose. 

We instead use the traditional gauge-fixed wall-source 
propagator and then compute either two-point functions 
for the J/ip and <fi states or four-point functions of the 
J/ip-cf) system using the wall-source operators. As we 
will show later, the two-hadron operator constructed by 
the wall-source single-hadron operators is automatically 
projected onto the trivial Ai irrep of point groups C nv 
in the two-hadron system. 



(> 



P a (fm) L 3 XT 


~ La (fm) 


csw 


Number of configs. 


1.9 0.0907(13) 32 3 x 64 


2.9 


1.715 


198 


TABLE II: Parameters of 2 + 1 flavor PACS-CS gauge 


configurations, generated 


using the Iwasaki 


gauge action and Wilson 


clover fermions, at = 156(7) MeV and vtik = 554(2) 









Flavor k v r s cb ce My (GeV) Reference 



Strange 


0.13640 


1.0 


1.0 


1.715 


1.715 


1.0749(33) 


[37] 


Charm 


0.10819 


1.2153 


1.2131 


2.0268 


1.7911 


3.0919(10) 





TABLE III: Parameters of clover fermions (strange) and the RHQ action (charm) used in this work and results of the vector 
meson masses. The strange quark mass is slightly off the physical point |37| . 



In the wall-source approach, the quark operator is 
summed over all spatial sites at the source time slice, 
where Coulomb gauge fixing is performed. This can 
be interpreted as zero-momentum-state projection in the 
quark-level kinematics. Under the twisted boundary con- 
ditions, however, the wall source applied to the Dirac ma- 
trix inversion with the hopping term defined in Eq. ([6| 
can be simply interpreted as the Fourier transformation 
to a momentum representation as 

^y(f^src) =^qe{x,t src )e- l(, - s/L -> q{p,t mc ) (17) 

X X 

with momentum given by p — 9/L. Thus, in this context, 
we directly construct the vector meson interpolating op- 
erator in three-dimensional momentum space from two 
wall sources with different twist angles 9 qf ^ 6 q , as 

0™(p,t src ) =Y,tf f ($,t STC )^Y,q' f ,(y,t SIC ) (18) 

x y 

with p = (9 qj — 6g f ,)/L. The superscript W stands for 
wall. 

Of course, when the same twisted boundary conditions 
such that 9 qj = 9q f , are imposed on both quark and an- 
tiquark operators, the resulting meson operator does not 
receive finite momentum. In the quarkonium case, such 
as J/tp and <fi, which are flavor-neutral mesons (/ = /'), 
9 qf = 9q f is a natural choice. Therefore, for the quarko- 
nium states, the twisted boundary condition trick does 
not seems to be applicable. However, we can practically 
choose 9 qf 7^ 9q f as partially twisted boundary conditions 
to construct the meson interpolating operator [33] . A 
price to pay is that, by construction, a disconnected di- 
agram contribution in quarkonium two-point functions 
is inevitably spoiled. In this context, the quarkonium 
states, which are flavor-neutral mesons, are fictitiously 
treated within this particular trick as if they were flavor- 
ful mesons consisting of two different flavors of the quark 
and antiquark with the same mass. 



In reality, the contributions from the disconnected di- 
agram in the vector channel are known to be negligibly 
small for strange and charm quarks in numerical sim- 
ulations [321 H3] in accordance with the mechanism of 
Okubo-Zweig-Iizuka suppression. Then, even if the case 
of the usual periodic boundary condition, namely, 9 = 0, 
is considered, the disconnected diagrams are often ne- 
glected in the J / ip and <f> meson spectroscopy since their 
evaluations are computationally expensive. Therefore, in 
this study, any disconnected diagram contribution (such 
as self-annihilations of both J/tp and <p) is not included in 
evaluating the two-point function of J / tp and <p mesons or 
the four-point functions of the J /tp-cp system in our sim- 
ulations regardless of whether we use the 9 = or 9 ^ 
case. For the twist angle, we choose \9 q \ ^ and \6q\ = 
or vice versa for the interpolating operator construction. 

At the sink we use a local bilinear operator as 

0%(x,t) = (f f (x,t)^q , f ,(x > t), (19) 

where L stands for local and / = /' should be taken 
for the J/tp and (f> states. As described in Eq. ([T]), un- 
der our choice of the partially twisted boundary condi- 
tions, a simple summation over the spatial sites on the 
sink operator automatically corresponds to the finite- 
momentum projection in terms of the Fourier transfor- 
mation as 0^(p,t) = ^2gO^(x,t) with either p = 9 q /L 

or —9q/L. We then construct two-point functions for a 
single quarkonium state with finite momentum p using 
the wall-source and local sink operators: 

G h (p, t, t SIC ) = \ J2&' H & WV^i-P, W 1 ) , (20) 
6 <=i 

where the superscript h stands for either the J/ip or <p 
states. We take an average over the spatial Lorentz in- 
dices i, which are related to the polarization direction of 
the states so as to obtain a possible reduction of statis- 
tical errors. Hereafter, we drop the label of both W and 
L on the single-hadron operators. 
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Oh 


-I Civ 


4- C2„ 


-I- C3,, 


p(s=0) 




Ax 


Ai 


Ai 


p(s=l) 


Tig 


A 2 © E 


a 2 © Si e s 2 


A 2 ® E 


p(s=2) 


Eg © T 2g 


(Ai e b x ) © (b 2 © e) 


(Ai © B2) © (Ai © A2 © Si) 


E © (Ai © £) 



TABLE IV: Subduction of the irreps of Oh, which are obtained from the decomposition of the direct product representation 
Tin <S> Ti u as the total spin of two vector particles, onto the little groups G<t v , C 2v and C$ v . 



Twist angle vector 


(0,0,0) 


(0,0, 9) 


(6,6,0) 


(0,9,9) 


(0,0, tt) 


(tt.it.O) 


(tt,tv,ty) 


Symmetry 


o h 




c 2v 


C3,, 




P>2h 


D 3d 


p(!=0) 


Alg 


Ax 


Ax 


Ax 


Alg 


A g 


Alg 


p(!=l) 




Ax (BE 


Ax © Bx © P 2 


Ax ®E 


A 2u © 


Bxu © -B2U © Bsu 


A 2u © E u 





TABLE V: Classification of the decomposition of the representations F*-' °- ) and 1 ' under each point group, which is 
associated with the choice of the twist angle vector 6. Point groups Oh, Dih, D 2 h and D^d have parity symmetry so that the 
trivial irrep Ai g (A 9 ) does not contain the I = 1 (P-wave) contribution. On the other hand, in the case of C^, v , C 2v and C 3v , 
where the inversion center is lost, the trivial irrep Ai contains both the I = (S-wave) and I = 1 (P-wave) contributions. 



Next we consider two-hadron interpolating operators 
for the J/ip-<f> system. In order to consider the CM sys- 
tem, we assign momentum p to the J/ip operator and the 
opposite momentum — p to the <p operator through the 
partially twisted boundary conditions imposed on both 
the charm quark and the strange quark having the same 
twist angle but with the opposite sign. The J/tp-f/) inter- 
polating operator for states with the relative momenta p 
of the J I system are constructed from the product of 
the single J/ip and (f> operators with opposite momentum: 



Q ij (p,t) = 0^(p,t)Of(-P,t + l), 



(21) 



where the subscripts i and j are spatial Lorcntz indices. 
To avoid the Fierz rearrangement of the J /ip-4> operator, 
the J/ip and <p operators are displaced by one time slice. 
Then, the four-point function of the J/ip-(p system can 
be constructed from the above defined operators: 



consider the partial wave (I) and spin (s) combination of 
the J/ip~4> system. To identify the J/ip-<p state labeled 
by the angular momentum I and spin s in continuum, we 
have to identify the corresponding irreducible represen- 
tations of the discrete rotation group on the lattice. 

The J/ip-(p system has three different spin states: spin- 
0, spin-1 and spin-2 states. However, the (2s + 1)- 
dimensional representation is, in general, reducible. 
In the case of 6 = 0, namely, p = 0, the spin-1 vector 
particle at rest is assigned to the T\ u irrep of the cubic 
group Oh [13 SS] ■ The total spin of two vector particles 
is given by the direct product representation T\ u © T\ u , 
which is decomposed into A ls © E g © T lg © T 2g , where 
Aig and E g j47] are one-dimensional and two-dimensional 
irreps respectively, while both T\ g and T 2g are three- 



dimensional irreps |45[ 146] . The following decomposition 
rules are then obtained: 



:) = (Qv(p,t)Qkl(~P,t SIc y). (22) 



Here we assume that the J / ip-(p system can be treated 
as nonrelativistic. The two-hadron operator defined in 
Eq. (21) can be described as the direct product of spin 



and momentum-space parts as 

Qi](p,t) = Sij © Q(p,t) 



(23) 



within the nonrelativistic approximation, where the 
partial-wave contributions are built in the spin- 
independent part, Q(p, t). This treatment would be justi- 
fied at least in the low-energy scattering region. We then 



r (s=o) 

p( S =l) 

p(s=2) 



= A 



= T, 



la- 



(24) 



T- 



2s 



E„. 



For the cases of spin and spin 1, there is a simple re- 
lationship between the cubic group representations and 
the continuum counterparts [481 . Both spin-0 and spin-1 
projections defined in the continuum theory remain valid 
for the corresponding irrep of the cubic group. We then 
obtain the respective four-point correlations for each ir- 
reducible representation of Oh as 



ljr A lg v 6 ' '■sic; 



r tJ/i>-4>/ + x \ 

T 2 g \ ' SrC ^ 



i=l j=l 

3 3 



j=l 3 = 1 

^EE{C(^M,g 



As one can easily check, a linear combination of two four- 
point correlation functions, which are projected onto the 
T 2g and E g irreps, reproduces the spin-2 projection in the 
continuum theory [35]. For the case of 9 = 0, where the 
cubic symmetry is satisfied, the spin-0, spin-1 and spin-2 
parts are orthogonal to each other. We then can easily 
determine the spin dependence of the J/tp-cj) interaction 
through four types of four-point correlation functions de- 
fined in Eqs. d25|-p8l). 



As we discussed in the previous section, for the cases 
of 9 0, the symmetry of the reciprocal lattice becomes 
C nv , which is the subgroup symmetry of the cubic group 
Oh, under the twisted boundary conditions. The decom- 



position rules ( 24 ) can be deduced by using the subduc- 
tion of the irrep of Oh to C± v , C 2v and C 3v as summa- 



rized in Table |IV| As can be seen in the new decompo- 
sition rules, the spin-0 operator as the trivial irrep A\ 
is inevitably mixed with the spin-2 operator under the 
twisted boundary conditions with 9^0. 

To clarify this point, let us consider the case of the 
point group Ci v , where symmetry is realized in the case 
of a twist angle vector of 9 = (0,0,6*). In this case, we 
may choose the corresponding spin-0, spin-1 and spin-2 
operators 

Qjy$-*®*) = Qu&t) + Q22@,t) + Qaa(p,t), 
Qj/U&V = Qi2(p,t)-Q 2 i(p,t), 

QjjU&t) = Qll(P,t) + Q22ip,t)-2Q 33 {p,t), 

which are guided by the knowledge in the continuum the- 
ory The spin-1 operator transforms according to the 
A 2 irrep of C^ v , while both spin-0 and spin-2 operators 
transform according to the Ai irrep. Although the spin-1 
operator is certainly orthogonal to the spin-0 and spin- 
2 operators, the latter ones are not orthogonal to each 
other in this case. Therefore, in practice, one may con- 
struct the 2x2 correlation matrix from Qjl^^ip, t) and 

Q^jJ^,_^(p,t) and then perform the variational method to 
disentangle the spin-0 and spin-2 contributions. 

Further reducing the symmetry from C^ v to C 2v , both 
the spin-1 and spin-2 parts are involved in the same ir- 



G 



{p = 0,t,i src )j , 



G%&+&= 0,t,i src ) - G^(P= 0,t,t SIC ) \ , 



2G^(p = 0, Msrc )}- 



(25) 
(26) 
(27) 
(28) 



r 



rep (A 2 ). Therefore, in this case, the disentanglement 
between the spin-1 and spin-2 contributions is also re- 
quired. Although the energy shift measured in the lower 
spin channel is slightly larger than the higher spin chan- 
nel, we do not find the appreciable spin dependence in the 
J /ip-(j) system through the calculation with 9 — 0, where 
the spin-0, spin-1 and spin-2 contributions are clearly dis- 
entangled as discussed previously. A similar finding in a 
quenched study of the J / ip-p interaction with the usual 
periodic boundary conditions was reported in Ref. [48] . 
This finding indicates that the spin-independent part of 
the J/ifj-tj) interaction dominates at least at low energies. 
Therefore, in this study, we will not resolve each spin 
contribution in simulations under the twisted boundary 
conditions, where the spin projection is somewhat com- 
plicated. Rather, we would like to resolve the partial- 
wave mixing between even-/ and odd-/ waves due to the 
twisted boundary conditions. 

Our aim in this study is to demonstrate the feasibility 
of our proposed approach, where both the .5- wave and P- 
wave phase shifts are extracted through the generalized 
Liischer finite size formula (14 1, as we will explain later. 
We focus on the spin-independent part of the J/tp~4> sys- 
tem for this purpose and thus use the spherical averaged 
(sav) J/tp-(f> operator for all types of the twisted angle 
vector 9: 



1 3 



(29) 



which transforms with respect to the rotation of the spin 
direction according to the trivial irreducible representa- 
tion of any point groups considered here. Then, the re- 
sulting four-point function turns out to be identical to 



the one defined in Eq. (25) for the irrep A\ g of Oh- 

Next we discuss the partial-wave contributions on the 
spherical averaged J/ip-<f> operator. The twisted bound- 
ary conditions make the symmetry of the reciprocal lat- 
tice break down to a subgroup symmetry of the cubic 
group Oh- Therefore, the (21 + l)-dimensional represen- 
tation rw becomes reducible and should be decomposed 
into the irreducible representations of the group C nv . For 
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instance, in the case of 9 = (0,0,0), there are four one- 
dimensional irreducible representations (A±, A%, B\, B2) 
and one two-dimensional irreducible representation (E). 
For the angular momenta I < 2, the resulting decompo- 
sitions are given [30H3II] as 



r C»=o) = Au 

r( ;=1 ) = Ai®e, 

r (?=2) = a 1 ®B 1 



(30) 



•Bi®E. 



This clearly indicates that the mixing between even-i and 
odd-Z wave contributions is not prohibited in some irreps. 
In this paper, we take advantage of this mixing in order 
to extract both S'-wave and P-wave phase shifts simul- 
taneously from the trivial A\ irrep of the two-hadron 
operator, which certainly contains both I = and / = 1 
contributions. There is nothing to change for either C^ v 
or where the trivial irrep A\ contains both S'-wave 
and P-wave contributions, like in C± v , as summarized in 
Table [V] As for the higher partial-wave contributions, we 
simply assume that they can be neglected. 

Here we note that the decomposition for the partial 
wave is almost similar to the case of the spin 
found in Table |IV[ except that the Ai irrep, instead of 
the A 2 irrep, appears in the second line. The rules are de- 
termined by the subduction of the irrep of the Oh group 
to the G± v group. The angular momentum I — 1 is as- 
signed to a vectorlike irreducible representation of Oh, 
which should be T\ u . The T\ u irrep of Oh is subdued 
to the A\ irrep of C& v 46 that appears in the second 
line of Eq. (30), while the subduction of the T\ g irrep is 



considered in Table |IV| for the total spin of two vector 
particles. 

In general, the operator that is projected onto some 
irrep ri rrop of the group G is simply given by 

9G R£G 



qa c — ' p 



REG 

where go is the number of elements R of the group G, 
while Xr iTrcp (R) is a character of R £ G for the irrep 
ri rre p. Here we consider that G is given by the groups 
C nv . In this case, the momentum p is invariant by any 
transformation R G G, namely, Rp — p. Therefore, we 
get 




(R)\Q(p,t), (31) 



Pr ittB „Q(p,t) 



where go = J2reg 1- Consequently, we observe that 



31) vanishes due to 



the right-hand side (r.h.s.) of Eq. 
T,ReGXr irlop ( R ) = 0, except for tl 
ial irrep A 1 , where ^2 ReG ^Ai( R ) 
fled. This means that the operator Q(p, t) is already 
irreducible and should transform according to the trivial 



re case of the triv- 
.9g 7^ is satis- 



irrep of the groups C nv [491 . This is a consequence of 
the wall-source propagator, where the quark operators 
are summed over all spatial sites at the source time slice. 
Thus, we can simply use the finite size formula of the Ai 
sector given in Eq. ( 14 ) since the two-hadron operator 



used in our simulations is constructed by the wall-source 
single-hadron operators. 



IV. NUMERICAL RESULTS 

We can read off the single-hadron energy from the two- 
point functions under the partial twisted boundary con- 
ditions: 



G h (p,t,t src ) 



-> e 



-E h (t-t st0 ) 



f»tsrc 



(32) 



where Eh = \Jp 2 + , with the rest mass Mh, which 
can be determined with the ordinary periodic boundary 
conditions (9 = 0). It is worth mentioning that Dirich- 
let boundary conditions are imposed for all quarks in 
the time direction in order to avoid wrapround effects, 
which are very cumbersome in systems of more than two 
hadrons. Each hadron mass is then obtained by fitting 
the corresponding two-point correlation function with a 
single exponential form. The (j) and J ftp masses are tab- 
ulated in Table EU 

Here we assume that each state (h = J/ip and 4>) on the 
lattice has the relativistic continuum dispersion relation, 
which is indeed enough to describe our data as shown in 
Fig. [2j On the other hand, the scattering momentum k 
of the J/if)-(f) two-particle system is defined through 



(33) 



where W denotes the total energy of the J/ip-cj) system in 
the CM frame. The CM energy can be determined from 
the large-< behavior of four-point functions 



Gj/^-^ip, t, t SIC ) 



t»ts 



-> e 



-W(t-W) 



(34) 



Next let us define the two-particle energy E measured 
from the threshold as 



E = W-(M J/ij +M <j> ). 



(35) 



The scattering momentum then can be represented in 
terms of E, 



_ y/E(E + 2M)(E + 2M. J/4 ,){E + 2Mj) 
k - 2{E + M) {3b) 

with the sum of the rest masses M = Mj/^ + M^ and the 
individual masses. Therefore, the precise determination 
of £ is a key point of the finite size analysis. Instead 
of trying to directly measure E, we first measure the 
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FIG. 2: Check of the dispersion relation for the (j> (upper 
panel) and J/ijj (lower panel) states. The vertical axis shows 
the momentum squared defined through the relativistic con- 
tinuum dispersion relation as pf on = E\ — M 2 for h = (f> and 
J/ip, while the horizontal axis is the momentum squared de- 
fined by the given twist angles as pf at = 8 2 / L 2 . By the linear 
fit (solid line) to data points calculated with various twist an- 
gles, the effective speed of light is obtained as c 2 ^ = 0.990(44) 
for the strange sector and c 2 a — 1.057(19) for the charm sec- 
tor. For comparison, the continuum dispersion relation is 
denoted as the dotted line in each panel. 



energy shift SE from the ratio of four-point and two-point 
functions with the partially twisted boundary conditions 



R 



G 



J/ip-<t> 



(p, t, t s 



t»t 3 



G J ^{p,t,t SIC )G<^{p,t + l,t aIC + l) 

(37) 



exp{~5E(t-t SJ . c )}, 



which reduces the statistical fluctuation due to a strong 
correlation between the denominator and numerator in 
the ratio. We then obtain the energy E from the follow- 
ing relation: 



E = SE + e 



J/i> 



(38) 



where = £7, — for h = J /if) and <p, which can 
be evaluated with the individual masses through the 
continuum-like dispersion relation. This procedure for 
the determination of E is similar to what was proposed 
inRef. Ell- 
in Fig. [3j we plot the effective energy shift defined by 



SE cS (p,t) = In 



R 



Rj/^-</,(p,t + i)'' 



(39) 



which should show a plateau for large Euclidean time 
(t t src ), for the case of a twist angle 9 — (0, 0, 9) with 
9 = 0.72 (rad) as a typical example. The quoted errors 
are estimated by a single elimination jackknifc method. 
An important observation is that the negative energy 
shift found in Fig. [3] implies the presence of an attrac- 
tive interaction between the J/ip and 4> states. 



FIG. 3: The effective energy shift SE in lattice units as a 
function of the time slice t, as in the case of a twist angle 
9 = (0, 0, 6) with 9 = 0.72 (rad). 



We find appropriate temporal windows, where the 
ground state dominance is satisfied, for fitting the ra- 
tio function defined in Eq. (37). Three horizontal solid 



lines represent the fit result with its 1 standard devia- 
tion obtained by a covariant single exponential fit over 
the range of 14 < t/a < 25, where two-point functions of 
the individual hadrons are also marginally dominated by 
the ground state of each hadron. We measure the energy 
shifts 8E for several twist angles and then observe that 
SE measured in the J/ip-cj) system is almost unchanged 
with the variation of the CM energy W. 

Once we obtain the scattering momentum k of the two- 
particle system, we can calculate the scattering phase 
shifts through a set of three finite size formulas ( 14 ) 



with three types of partially twisted boundary conditions, 
namely, [001], [110] and [111]. These formulas contain 
two unknowns, namely the S- wave and P-wave phase 
shifts. In order to obtain both phase shifts at various 
scattering momenta, we propose the following strategy. 

First, we calculate the 5-wave phase shift So(k) in the 
case of four special twist angles, 9 = (0,0,0), (0,0, n), 
(7r,7r, 0) and (7r,7r,7r), where the parity symmetry re- 
mains preserved. Since there is no unwanted mixing 
between even-Z and odd-/ wave contributions, the finite 
size formula reduces to either the Luscher formula or 
the Rummukainen-Gottlieb— type formula without the 
Lorentz factor (since 7 = 1). In Fig. [4j we plot results 
of fccotJrj, which are obtained from the data calculated 
with specific twist angles, 9 = (0,0,0), (0,0, it), (tt,tt, 0) 
and (it, n, 7t) as a function of k 2 in lattice units. We ob- 
serve that /ccot(5o monotonically increases as k 2 increases 
and there is no nonanalytic behavior in the range cov- 
ered here. Therefore, we simply interpolate the four data 
points by the quadratic function of k 2 using the jackknife 
method. 

Next, using the above information of the S- wave phase 
shift <5q, we evaluate cotc5i from the data of the scattering 
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momentum k 2 calculated with a twist angle vector, 
(6*, 9, 9) through the formula 



cotdo(fc) - M l ss '(q) 



woo(q) + zVGlm {w 2 2(q)} 



9w 2 (q) 



cot S (k) - woo(q) 
(40) 



The results obtained for the P-wave phase shift S\ are 
plotted as fc 3 cot(5i vs k 2 in Fig. 5. Although a slight 
but monotonic increase is observed, the results of fc 3 cot<5i 
show only very mild k 2 dependence. This observation 
suggests that the effective-range expansion is valid in this 
k 2 region. 




0.005 0.01 0.015 0.02 0.025 0.03 

k 2 



FIG. 4: fccot^o as a function of k 2 in lattice units. Open 
circles are calculated with twist angles, 9 = (0,0,0), (0,0, it), 
(n,n,0) and (n,n,n). A solid curve represents an interpola- 
tion of four data points using the effective range expansion 
up to the order of fc 4 . 

In this work, we use five different angles, 9 = 1.84, 1.92, 
2.14, 2.35 and 2.56 (rad), to choose k 2 points that range 
from (j) to2(|j) in lattice units. This is because when 
the value of 9 becomes very small, the resulting k 2 is close 
to the threshold, k 2 ~ 0. As we will discuss later, the P- 



wave mixing contributions to the finite size formula ( 14 ) 



become negligible in the region of k 2 < f~) so that the 
term of cot(5o — wqq is almost vanishing. This gives rise to 
numerical difficulties in this region for extracting Si from 
Eq. (40 ), where the second term on the right-hand side is 



very sensitive to how precisely the scattering momentum 
k 2 is determined. 

To avoid the above numerical problem, we evaluate 
fc 3 cot<5 1 up to k 2 ~ 0.008 from 0.02 in lattice units, where 
the term of cot<5o — woo is determined precisely enough, 
and we then extrapolate the value of fc 3 cot<5i to the re- 
gion k 2 ~ by the polynomial function of k 2 , which cor- 
responds to the effective-range expansion 50J . Thanks 
to very mild k 2 dependence, the linear function of k 2 is 




0.02 



FIG. 5: fc 3 cot5i as a function of k 2 in lattice units. Open 
diamonds are calculated with a twist angle vector 9 = (9, 9, 9) 
with five different angles 9. For evaluation of cot<5i, Eq. (401 
is used with an input of cot^o, which was determined in Fig.|4j 

enough to extrapolate the data of /c 3 cot<5 1 from the cur- 
rent set of points toward the low k 2 region, since the 
higher-order terms in the k 2 expansion should be sup- 
pressed near the threshold. 

We thus obtain one of the threshold parameters for the 
P wave, the scattering volume oi, which is defined as an 
inverse of fc 3 cot(5i determined at k = 0: 



ai 



tanJi 



fc 3 



0.0348 ± 0.0081 [fm 3 



(41) 



fc=0 



through the linear fit on the data of fc 3 cot<5i. The result 
of ai indicates that the P-wave interaction of the J /ip-<j) 
system is weakly attractive. 

Combined with this P-wave information, we also eval- 
uate fccot<5o at very low energies using the data calculated 
under other twisted boundary conditions with a twist an- 
gle vector 9 = (0,0,0), through the finite size formula 



COtSo(fe) = M [ gf ] {q) + 



w o{q) 



\Mf P %W 



cot(5i(fc)-A^ 1] ((7) 



3w( (q) 



cot 5i(k) - w 00 (q) - 2w 20 (q) ' 



(42) 



2.16 and 

if 



Here we use three different angles, 9 = 1.44 
2.88 (rad), to choose k 2 points that range from to 
in lattice units. 

Figure [6] shows the results of fccot<5o . We have carried 
out correlated \ 2 fits on all seven data displayed in Fig. [6] 
by using the form of the effective-range expansion: 



kcotSo{k) = 



«o 



r k 2 



N 
n=2 



k 2 



(43) 



up to N = 4, which corresponds to the order of fc 8 . We 
recall that the effect of Z?-wave contributions is, strictly 
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O (0,0,0), (0,0,11), (Jt.it.O), (7t,7t,Jt) 
□ (0,0,9) with 3 different angles 6 
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j. 
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FIG. 6: fccot^o as a function of k 2 in lattice units, with 
9 — (0, 0, 9) data. The same data that appeared in Fig. [4] 
are shown as open circles, while open squares stand for new 
data calculated with another twist angle vector 9 = (0, 0, 9) 
with three different twist angles. For the evaluation of cot Jo, 
Eq. (42 1 is used with an input of cot(5i , which was determined 
in Fig[5| A dashed curve with a band indicates the final fit 
result with 1 standard deviation, which is obtained from re- 
newed fitting on three new data together with four data. 



speaking, not negligible at the order of k 4 ; the coefficient 
V2 in the k 4 term, which corresponds to the scattering 
volume prj) , may suffer from systematic uncertainties due 
to such unknown effects of higher partial-wave mixing. 
However, we will quote the coefficient v 2 as a reference 
value of the scattering volume as well as the scattering 
length ao and the effective range ro, which are safely read 
off from the coefficients of the k° and k 2 terms later. 

The stability of the fit results has been tested against 
either the number of fitted data points or the number of 



polynomial terms for a given order N defined by Eq. ( 43 1 



The best fit is drawn to fit all seven data points using the 
fitting form ( 43 1 with N — 4 in Fig. [6j We then obtain 
the following results: 



Oo 

ro 
Pro 



0.242 ±0.041 [fm], 
4.712 ± 1.727 [fm], 
2.50 ±1.64 [fm 3 ], 



where the quoted errors represent only the statistical er- 
rors given by the jackknife analysis. The effective range 
ro receives a rather large error. This is because the size 
of the effective range is much larger than that of the 
scattering length, r 3> do. In this particular case, the 
coefficient of the linear term in momentum space is in- 
fluenced even by small statistical fluctuations at lower 
energy data points. More statistics are clearly needed to 
reduce the statistical error on the effective range. But, 
we would like to stress here that our approach shows fea- 
sibility for determining the effective range as well as the 
scattering length. 

Figure [7] shows S-wave (left panel) and P-wave (right 
panel) scattering phase shifts for the J/if>-<l> system. The 



phase shifts observed here are positive reflecting the at- 
tractive interaction between the J / ip and <f> states in both 
channels. In each panel, a dashed curve represents the 
fit result guided by the effective range expansion applied 
to all data points except for open left-triangle symbols, 
which are additionally calculated with the other twist 
angle vector 6 = (6,6,0). The fit result of the P-wave 
phase shift is used as prior information in the S'-wave cal- 
culation for new data, and vice versa. New data points 
in both panels show an excellent agreement with the fit 
curves. This means that the results obtained from our 
data analysis developed here pass the consistency test. 

Our results exhibit typical behaviors of weak attraction 
(small scattering length and small scattering volume) in 
S'-wave and P-wave phase shifts at low energies. Unfor- 
tunately, there is no resonance structure associated with 
the F(4140) resonance against what we expected. Al- 
though the CDF experiment has reported evidence for 
a new narrow resonance, F(4140) around a few 10 MeV 
above the J/ip-'P threshold, it seems that our results are 
consistent with null results in the Belle and LHCb exper- 
iments. 

Next, we give some remarks on how large the effect of 
the higher partial-wave mixing remains at low energies. 
Although the lower partial wave becomes dominant at 
low energies due to Si(p) oc aip 2l+1 , the P-wave mixing 
induced by the usage of the partially twisted boundary 
conditions should be suppressed in the limit of 6 —> or 
L — >• 00. If we assume |Af^p| ~ at low energies, the 
finite size formula then reduces to the original Luscher's 
formula 

coU (fc) = w 00 (q) = -^Z^il-q 2 ) (44) 
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with the twist angle vector 6. 

Figure [8] exposes the size of the systematic error in the 
analysis where the P-wave mixing is ignored. Open sym- 
bols are calculated through the approximated formula 
(J44[ ) for data taken with two types of twist angle vectors, 
6 = (0,0,61) and {6,6,6). In Fig. § a solid curve corre- 
sponds to the S-wave phase shift in the full analysis as 
described previously, while a dashed curve represents the 
P-wave phase shift in the same analysis. When the P- 
wave phase shift is smaller enough than the S-wave phase 
shift, S ^> Si, the difference between the results from the 
full and approximated analyses gets smaller, especially 
near the threshold, as we expected. However, the validity 
region of the approximation, such that |A4gp| = 0, is lim- 
ited in the vicinity of the threshold even for such weakly 
interacting systems, where both the scattering length and 
the scattering volume are observed to be small. 

From this observation, our ignorance of the higher par- 
tial wave (I > 2) in the full analysis would be acceptable, 
at least up to a few 10 MeV above the threshold, where 
the P-wave mixing is safely ignored. We therefore ex- 
pect that our null result for the F(4140) resonance in 
this vicinity in both S- and P-wave channels remains un- 
changed even if the D-wsve mixing is taken into account 
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FIG. 7: S'-wave (left panel) and P-wave (right panel) scattering phase shifts as a function of the two-particle energy E, which 
is measured from the J/ip-(f> threshold. In each panel, a dashed curve represents the fit result guided by the effective range 
expansion applied to all data points except for open left-triangle symbols, which are additionally calculated with the other 
twist angle vector 8 = (8,8,0). 



in the finite size formula. Nevertheless, we plan to de- 
velop our analysis, with the finite size formula extended 
to treat the higher- wave mixing up to I = 2, in a sepa- 
rate publication. As such, the formula has already been 
considered for the case of general moving-frame compu- 
tation 133 . 



Finally, we make a comment about the usage of the 
partially twisted boundary conditions, which are imposed 
only on valence quarks and then may induce some sys- 
tematic error due to the difference of boundary conditions 
between dynamical and valence quarks. This is only an 
issue for the strange quark in the J/ip-cp system, since 
the charm quark is treated as a quenched approxima- 
tion. The authors of Ref. [25] show that the finite volume 
correction due to the partially twisted boundary condi- 
tions is exponentially suppressed as the spatial extent L 
increases. The coefficient of exponential-type finite size 
corrections may be characterized by a mass of the lightest 
state mediated between two-particles. In the J/ip-fp sys- 
tem, a single meson exchange such as one pion exchange 
is forbidden, unlike a typical hadron-hadron interaction, 
since the J/ip and <f> states do not carry any isospin and 
also do not share the same quark flavor. Therefore, in- 
stead of the complex nature originating from quark ex- 
changes, multigluon exchanges become dominant in the 
J/ip-(j) system. An effective description for the interac- 
tion between the J/tp and <fi may be given by the soft 
Pomeron exchange, which carries the quantum number of 
the vacuum |51j . In this sense, the difference of bound- 
aries between dynamical quarks and valence quarks in 
this study would be negligibly small since the mass of 
the Pomeron is typically of the order of 500—600 MeV, 
which is, phenomenologically, a good agreement with 
what is observed through the charmonium-nucleon po- 
tential EHE21- 
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FIG. 8: Scattering phase shifts as a function of the two- 
particle energy E. Ope n sy mbols are calculated through the 
approximated formula (441, where the P-wave mixing is ig- 



nored, for data taken with two types of twist angle vectors, 
8 — (0, 0, 8) and (8, 8, 8). Solid and dashed curves correspond 
to the S'-wave and P-wave phase shifts given in the full anal- 
ysis, as described in Fig. [7] 



V. SUMMARY 

We have considered how to develop the Liischer finite 
size approach under the partially twisted boundary con- 
dition. The twisted boundary condition allows us to treat 
any small momentum on the lattice through the variation 
of the twist angle, continuously. On the other hand, when 
the twisted boundary condition is introduced, the cubic 
symmetry is broken down to some subgroup symmetry, 
where the center of inversion symmetry is unfortunately 
lost. Accordingly, even-l and odd-Z partial waves are in- 
evitably mixed together in the extended Liischer finite 
size formula. 
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We propose to take advantage of these unique prop- 
erties of the twisted boundary condition in order to ex- 
tract both the S- wave and P-wave scattering phase shifts 
from the two-particle energy with the total momentum 
\P\ = calculated in a single finite box with various 
types of partially twisted boundary conditions. We then 
demonstrate the feasibility of a new approach to exam- 
ine the existence of near-threshold and narrow resonance 
states, whose characteristics are observed in many of the 
newly discovered charmonium-like XYZ mesons. 

As an example, we choose low-energy J/tp-(j> scatter- 
ings to search for the Y(4140) resonance, which is ob- 
served around a few 10 MeV above the J/ip-(f> thresh- 
old by the CDF experiment at Fermilab through the 
B + — > J/ip<frK + decay process from pp collisions. Our 
simulations are performed in 2+1 flavor dynamical lat- 
tice QCD using the PACS-CS gauge configurations at the 
lightest pion mass, m, = 156 MeV, with a relativistic 
heavy-quark action for the charm quark. We successfully 
obtain the P-wave phase shift of the J/il>-<j> scattering as 
well as the S- wave phase shift near the threshold. 

Our results exhibit typical behaviors of weak attrac- 
tion (small scattering length and small scattering vol- 
ume) in both .S'-wave and P-wavc channels at low en- 
ergies. There is no resonance structure associated with 
the Y"(4140) resonance, contrary to what we expected. 
Therefore, our results are consistent with null results 



in the Belle and LHCb experiments. Instead of find- 
ing near-threshold and narrow resonance states, we suc- 
ceed in extracting model-independent information of the 
low-energy J/t/j-<f> interaction, such as the scattering 
length ao = 0.242 ± 0.041 fm and the effective range 
ro = 4.71 ± 1.73 fm for the 5-wave and the scattering 
volume ai = 0.035 ± 0.008 fm 3 for the P-wave, within 
our new approach. We plan to apply this new method 
to nucleon-nucleon scatterings and also D-K and D*-K 
scatterings. The former is a key target in the original 
proposal of twisted boundary conditions [12] , while the 
latter would give some insight into the structure of the 
D s o(2317) and P s i(2460) resonances. 



Acknowledgments 

We would like to thank T. Hatsuda and T. Kawanai 
for fruitful discussions. This work is supported by MEXT 
Grant-in- Aid for Scientific Research on Innovative Areas 
(No.2004:20105003) and JSPS Grants-in-Aid for Scien- 
tific Research (C) (No. 23540284). Numerical calcula- 
tions reported here were carried out on the T2K super- 
computer at ITC, University of Tokyo and also at CCS, 
University of Tsukuba. 



[1] M. Liischer, Commun. Math. Phys. 105, 153 (1986). 

[2] M. Liischer, Nucl. Phys. B 354, 531 (1991). 

[3] S. R. Beane, K. Orginos and M. J. Savage, Int. J. Mod. 

Phys. E 17, 1157 (2008) and references therein. 
[4] Signatures of S'-wave bound state formation in finite 

volume can also be discussed on the basis of Luscher's 

method 00. 

[5] S. R. Beane, P. F. Bedaque, A. Parreno and M. J. Savage, 

Phys. Lett. B 585, 106 (2004). 
[6] S. Sasaki and T. Yamazaki, Phys. Rev. D 74, 114507 

(2006) . 

[7] K. Rummukainen and S. A. Gottlieb, Nucl. Phys. B 450, 
397 (1995). 

[8] N. H. Christ, C. Kim and T. Yamazaki, Phys. Rev. D 

72, 114506 (2005). 
[9] C. h. Kim, C. T. Sachrajda and S. R. Sharpe, Nucl. Phys. 

B 727, 218 (2005). 
[10] X. Feng, X. Li and C. Liu, Phys. Rev. D 70, 014505 

(2004). 

[11] X. Li et al. [CLQCD Collaboration], JHEP 0706, 053 

(2007) . 

[12] P. F. Bedaque, Phys. Lett. B 593, 82 (2004). 

[13] R. G. Newton, "Scattering Theory of Waves and Parti- 
cles", 2nd ed. (Springer, New York, 1982). 

[14] S. Ozaki and S. Sasaki, PoS LATTICE2012, 160 (2012). 

[15] N. Brambilla, S. Eidelman, B. K. Heltsley, R. Vogt, 
G. T. Bodwin, E. Eichten, A. D. Frawley and A. B. Meyer 
et al, Eur. Phys. J. C 71, 1534 (2011). 

[16] S. Godfrey and S. L. Olsen, Ann. Rev. Nucl. Part. Sci. 
58, 51 (2008). 



[17] A. Bondar et al. [Belle Collaboration], Phys. Rev. Lett. 

108, 122001 (2012). 
[18] I. Adachi [Belle Collaboration] , |arXiv: 1 105.4583| [hep-ex] . 
[19] K. F. Chen et al. [Belle Collaboration], Phys. Rev. Lett. 

100, 112001 (2008). 
[20] K. -F. Chen et al. [Belle Collaboration], Phys. Rev. D 

82, 091106 (2010). 
[21] T. Aaltonen et al. [CDF Collaboration], Phys. Rev. Lett. 

102, 242002 (2009). 
[22] T. Aaltonen et al. [The CDF Collaboration], 

|arXiv:1101.6058| [hep-ex]. 
[23] C. P. Shcn et al. [Belle Collaboration], Phys. Rev. Lett. 

104, 112004 (2010). 
[24] R. Aaij et al. [LHCb Collaboration], Phys. Rev. D 85, 

091103 (2012). 

[25] C. T. Sachrajda and G. Villadoro, Phys. Lett. B 609, 73 
(2005). 

[26] G. M. de Divitiis, R. Petronzio and N. Tantalo, Phys. 

Lett. B 595, 408 (2004). 
[27] P. A. Boyle et al., JHEP 0807, 112 (2008). 
[28] C. H. Kim and C. T. Sachrajda, Phys. Rev. D 81, 114506 

(2010) . 

[29] Z. Davoudi and M. J. Savage, Phys. Rev. D 84, 114502 

(2011) . 

[30] Z. Fu, Phys. Rev. D 85, 014506 (2012). 

[31] L. Leskovec and S. Prelovsek, Phys. Rev. D 85, 114507 

(2012) . 

[32] M. Doring, U. G. Meissner, E. Oset and A. Rusetsky, 

Eur. Phys. J. A 48, 114 (2012). 
[33] M. M. Gockeler, R. Horsley, M. Lage, U. -G. Meiss- 



15 



ner, P. E. L. Rakow, A. Rusetsky, G. Schierholz and 
J. M. Zanotti, Phys. Rev. D 86, 094513 (2012). 

[34] T. Ka wanai and S. Sasaki, PoS LATTICE 2010, 156 
(2010) [arXiv: 101 1 .1322] [hep-lat] ] . 

[35] The Lorentz transformation is characterized by the ve- 
locity v = P/E, where E and P are the total energy and 
the total momentum, respectively, and also the Lorentz 
boost factor 7 = l/yl — W. 

[36] T. Yamazaki et al. [CP-PACS Collaboration], Phys. Rev. 
D 70, 074513 (2004). 

[37] S. Aoki et al. [PACS-CS Collaboration], Phys. Rev. D 
79, 034503 (2009). 

[38] http://www.usqcd.org/ildg (ILDG) and 

http://www.jldg.org (JLDG). 

[39] A. X. El-Khadra et al, Phys. Rev. D 55, 3933 (1997). 

[40] S. Aoki, Y. Kuramashi and S. I. Tominaga, Prog. Theor. 
Phys. 109, 383 (2003). 

[41] Y. Kayaba et al. [CP-PACS Collaboration], JHEP 0702, 
019 (2007). 

[42] T. Kawanai and S. Sasaki, Phys. Rev. D 85, 091503 
(2012). 

[43] C. McNeile et al. [UKQCD Collaboration], Phys. Rev. D 
70, 034506 (2004). 



[44] P. de Forcrand et al. [QCD-TARO Collaboration], JHEP 
0408, 004 (2004). 

[45] D. C. Moore and G. T. Fleming, Phys. Rev. D 73, 014504 
(2006) [Erratum-ibid. D 74, 079905 (2006)]. 

[46] J. Foley, J. Bulava, Y. -C. Jhang, K. J. Juge, D. Lenkner, 
C. Morning star and C. H. W ong, PoS LATTICE 2011, 
120 (2011) [arXiv: 1205.4223 [hep-lat]]. 

[47] The index g (gerade) or u (ungerade) denotes the exis- 
tence of an inversion center and indicates even or odd 
parity. 

[48] K. Yokokawa, S. Sasaki, T. Hatsuda and A. Hayashigaki, 

Phys. Rev. D 74, 034504 (2006). 
[49] For the case of 9 — (0, 0, 0), where G is given by Oh, the 

condition of Rp = p is automatically satisfied since p — 0. 
[50] The quadratic coefficient in the k 4 term should be un- 

physical, since the effect of higher partial-wave (D-wave) 

contributions is not negligible at the order of fe 4 . 
[51] P. V. Landshoff and O. Nachtmann, Z. Phys. C 35, 405 

(1987). 

[52] T. Kawanai and S. Sasaki, Phys. Rev. D 82, 091501 
(2010). 



